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Figure 1: Partial functional correspondence between two pairs of shapes with large missing parts. For each pair we show the 
matrix C representing the functional map in the spectral domain, and the action of the map by transferring colors from one 
shape to the other. The special slanted-diagonal structure of C induced by the partiality transformation is first estimated from 
spectral properties of the two shapes, and then exploited to drive the matching process. 


Abstract 

In this paper, we propose a method for computing partial functional correspondence between non-rigid shapes. We 
use perturbation analysis to show how removal of shape parts changes the Laplace-Beltrami eigenfunctions, and 
exploit it as a prior on the spectral representation of the correspondence. Corresponding parts are optimization 
variables in our problem and are used to weight the functional correspondence; we are looking for the largest 
and most regular (in the Mumford-Shah sense) parts that minimize correspondence distortion. We show that our 
approach can cope with very challenging correspondence settings. 

Categories and Subject Descriptors (according to ACM CCS): 1.3.5 [Computer Graphics]: Computational Geometry 
and Object Modeling—Shape Analysis 


1. Introduction 

The problem of shape correspondence is one of the most 
fundamental problems in computer graphics and geometry 
processing, with a plethora of applications ranging from 
texture mapping to animation [BBK06, KLCF10, KLF11, 
VKZHCOll]. A particularly challenging setting is that of 
non-rigid correspondence , where the shapes in question are 
allowed to undergo deformations, which are typically as¬ 
sumed to be approximately isometric (such a model appears 
to be good for, e.g. , human body poses). Even more chal¬ 
lenging is partial correspondence , where one is shown only 
a subset of the shape and has to match it to a deformed full 
version thereof. Partial correspondence problems arise in nu¬ 
merous applications that involve real data acquisition by 3D 
sensors, which inevitably lead to missing parts due to occlu¬ 
sions or partial view. 


Related work. For rigid partial correspondence problems, 
arising e.g., in 3D scan completion applications, many ver¬ 
sions of regularized iterative closest point (ICP) approaches 
exist, see for example [AMCO08, ART 15]. Attempts to ex¬ 
tend these ideas to the non-rigid case in the form of non- 
rigid or piece-wise rigid ICP have been explored in recent 
years [LSP08]. By nature of the ICP algorithm, these meth¬ 
ods rely on the assumption that the given shapes can be 
placed in approximate rigid alignment to initiate the match¬ 
ing process. As a result, they tend to work well under small 
deformations (e.g., when matching neighboring frames of a 
sequence), but performance deteriorates quickly when this 
assumption does not hold. 

For the non-rigid setting, several metric approaches cen¬ 
tered around the notion of minimum distortion correspon¬ 
dence [BBK06] have been proposed. Bronstein et al. [BB08, 
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BBBK09] combine metric distortion minimization with op¬ 
timization over matching parts, showing an algorithm that 
simultaneously seeks for a correspondence and maximizes 
the regularity of corresponding parts in the given shapes. 
Rodola et al. [RBA*12] subsequently relaxed the regular¬ 
ity requirement by allowing sparse correspondences, and 
later introduced a mechanism to explicitly control the de¬ 
gree of sparsity of the solution [RTH* 13]. Finally, in [SY14] 
the authors proposed a voting-based formulation to match 
shape extremities, which are assumed to be preserved by the 
partiality transformation. Being based on spectral features 
and metric preservation, the accuracy of the aforementioned 
methods suffers at high levels of partiality, where the compu¬ 
tation of these quantities becomes unreliable due to bound¬ 
ary effects and meshing artifacts. Furthermore, these meth¬ 
ods suffer from high computational complexity and gener¬ 
ally provide only a sparse correspondence. 

Pokrass et al. [PBB13] proposed a descriptor-based par¬ 
tial matching approach where the optimization over parts is 
done to maximize the matching of bags of local descriptors. 
The main drawback of this approach is that it only finds simi¬ 
lar parts, without providing a correspondence between them. 
Windheuser et al. [WSSC11] formulated the shape match¬ 
ing problem as one of seeking minimal surfaces in the prod¬ 
uct space of two given shapes; the formulation notably al¬ 
lows for a linear programming discretization and provides 
guaranteed continuous and orientation-preserving solutions. 
The method was shown to work well with partial shapes, but 
requires watertight surfaces as input ( e.g ., via hole filling). 
Brunton et al [BWW* 14] used alignment of tangent spaces 
for partial correspondence. In their method, a sparse set of 
correspondences is first computed by matching feature de¬ 
scriptors; the matches are then propagated in an isometric 
fashion so as to cover the largest possible regions on the two 
shapes. Since the quality of the final solution directly de¬ 
pends on the initial matches, the method is understood as 
a “densification” method to complement other sparse ap¬ 
proaches. Other recent works include the design of robust 
descriptors for partial matching [vKZH13]. In the context of 
collections of shapes, partial correspondence has been con¬ 
sidered in [YKTS* 11,CGH14,CRA* 16]. 

All the aforementioned works are based on the notion of 
point-wise correspondence between shapes. Recently, Ovs- 
janikov et al. [OBCS*12] proposed the functional maps 
framework, in which shape correspondence is modeled as 
a linear operator between spaces of functions on the shapes. 
The main advantage of functional maps is that finding cor¬ 
respondence boils down to a simple algebraic problem, 
as opposed to difficult combinatorial-type problems aris¬ 
ing in, e.g., the computation of minimum-distortion maps. 
While several recent works showed that functional maps 
can be made resilient to missing parts or incomplete data 
[HWG14, KBBV15], overall this framework is not suitable 
for dealing with partial correspondence. 


Contribution. In this paper, we propose an extension to the 
functional correspondence framework to allow dealing with 
partial correspondence. Specifically, we consider a scenario 
of matching a part of a deformed shape to some full model. 
Such scenarios are very common for instance in robotics 
applications, where one has to match an object acquired 
by means of a 3D scanner (and thus partially occluded) 
with a reference object known in advance. We use an ex¬ 
plicit part model over which optimization is performed as 
in [BB08,BBBK09], as well as a regularization on the spec¬ 
tral representation of the functional correspondence account¬ 
ing for a special structure of the Laplacian eigenfunctions as 
a result of part removal. Theoretical study of this behavior 
based on perturbation analysis of Laplacian matrices is an¬ 
other contribution of our work. We show experimentally that 
the proposed approach allows dealing with very challenging 
partial correspondence settings; further, we introduce a new 
benchmark to evaluate deformable partial correspondence 
methods, consisting of hundreds of shapes and ground-truth 
information. 

The rest of the paper is organized as follows. In Section 
2, we review the basic concepts in the spectral geometry and 
describe the functional correspondence approach. Section 3 
studies the behavior of Laplacian eigenfunctions in the case 
of missing parts, motivating the regularizations used in the 
subsequent sections. Section 4 introduces our partial corre¬ 
spondence model, and Section 5 describes its implementa¬ 
tion details. Section 6 presents experimental results, and fi¬ 
nally, Section 7 concludes the paper. 

2. Background 

In this paper, we model shapes as compact connected 2- 
manifolds M, possibly with boundary dM. Given f,g : 
M —> K some real scalar fields on the manifold, we de¬ 
fine the standard inner product ( f,g)M = Im f( x )s( x )d x > 
where integration is done using the area element in¬ 
duced by the Riemannian metric. We denote by L 2 (M) = 
{/ : M R | (/,/) Af < oo} the space of square-integrable 
functions on M . 

The intrinsic gradient V m! and the positive semi- 
definite Laplace-Beltrami operator A = — div^t (S M.f) 

generalize the notions of gradient and Laplacian to man¬ 
ifolds. The Laplace-Beltrami operator admits an eigen- 
decomposition 

4m W*) = * £ int {M) (1) 

(V.M<t>i(x) >»(*))= 0 xedM, (2) 

with homogeneous Neumann boundary conditions (2) if M 
has a boundary (here n denotes the normal vector to the 
boundary), where 0 = A-i < X 2 < ... are eigenvalues and 
(|)l, (f> 2 , • • • are the corresponding eigenfunctions (or eigen¬ 
vectors). The eigenfunctions form an orthonormal basis 
on L 2 (J W), i.e., (§i,§j)M = &ij> generalizing the classical 
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Fourier analysis: a function / £ L 2 (M) can be expanded 
into the Fourier series as 

m = o) 

*>i 

Functional correspondence. Let us be now given two man¬ 
ifolds, J\f and M. Ovsjanikov et al. [OBCS*12] proposed 
modeling functional correspondence between shapes as a 
linear operator T : L 2 (AT) —» L 2 (Ad). One can easily see that 
classical vertex-wise correspondence is a particular setting 
where T maps delta-functions to delta-functions. 

Assuming to be given two orthonormal bases {(])/}/>i and 
{V*'}i>l on L 2 (JV) and L 2 (Ad) respectively, the functional 
correspondence can be expressed w.r.t. to these bases as fol¬ 
lows: 

Tf = 

i> 1 i> 1 

= £ ( 4 ) 

<7>l v ---' 

C ij 

Thus, T amounts to a linear transformation of the Fourier 
coefficients of / from basis {())/}/>i to basis {\|//}/>i, which 
is captured by the coefficients Qy. Truncating the Fourier se¬ 
ries (4) at the first k coefficients, one obtains a rank-/; ap¬ 
proximation of T, represented in the bases {<])/,\|//}/>i as a 
kxk matrix C = (oy). 

In order to compute C, Ovsjanikov et al. [OBCS*12] 
assume to be given a set of q corresponding functions 
{/l, • • • ,/?} C L 2 (A/") and {gi,... ,#$} C L 2 (M). Denoting 
by a ij = (fj><h)u and hj = {gj,^i)M die k x q matrices 
of the respective Fourier coefficients, functional correspon¬ 
dence boils down to the linear system 

CA = B (5) 

If q > k, the system (5) is (over-)determined and is solved in 
the least squares sense to find C. 

Structure of C. We note that the coefficients C depend on 
the choice of the bases. In particular, it is convenient to use 
the eigenfunctions of the Laplace-Beltrami operators of J\f 
and Ad as the bases {(j) z , V|//};>i; truncating the series at the 
first k coefficients has the effect of ‘low-pass’ filtering thus 
producing smooth correspondences. In the following, this 
will be our tacit basis choice. 

Furthermore, note that the system (5) has qk equations and 
k 2 variables. However, in many situations the actual number 
of variables is significantly smaller, as C manifests a cer¬ 
tain structure which can be taken advantage of. In particular, 
if A f and Ad are isometric and have simple spectrum (i.e., 
the Laplace-Beltrami eigenvalues have no multiplicity), then 
T§i = =L\|/;, or in other words, Qy = ±5 y. In more realis¬ 
tic scenarios (approximately isometric shapes), the matrix C 
would manifest a funnel-shaped structure, with the majority 
of elements distant from the diagonal close to zero. 


Discretization. In the discrete setting, the manifold A f is 
sampled at n points x\,...,x n which are connected by edges 
E = E{ U£b and faces F, forming a manifold triangular mesh 
(V,E,F). We denote by E x and E^ the interior and bound¬ 
ary edges respectively. A function on the manifold is repre¬ 
sented by an ^-dimensional vector f = (/(*i),...,/(*«)) T . 
The discretization of the Laplacian takes the form of an n x n 
sparse matrix L = — S -1 W using the classical cotangent for¬ 
mula [Mac49, Duf59, PP93], 


(cot QLij + cotP ij)/2 ij £ E x ; 
(cotcqy)/2 ij CE h ; 

-LkfiWik i = j\ 

0 else; 


( 6 ) 


where S = diag(si,..., s n ), = \ ZjhijkeF s ijk denotes the 
local area element at vertex /, s^k denotes the area of triangle 
i jk, and oqy, P, y denote the angles Zikj, Z jhi of the triangles 
sharing the edge ij (see Fig. 2). 




Figure 2: Discretization of the Laplace-Beltrami operator 
on a triangular mesh for interior edges (green, left) and 
boundary edges (red, right). 


The first k eigenfunctions and eigenvalues of the Lapla¬ 
cian are computed by performing the generalized eigen- 
decomposition W<£ = SOA, where <£ = (4> l5 ... is an 
n X k matrix containing as columns the discretized eigen¬ 
functions and A = diag(?ti,..., Xjj) is the diagonal matrix of 
the corresponding eigenvalues. The computation of Fourier 
coefficients is performed by a = 4> T Sf. 


3. Laplacian eigenvectors and eigenvalues under 
partiality 

When one of the two shapes has missing parts, the assump¬ 
tion of approximate isometry does not hold anymore and a 
direct application of the method of Ovsjanikov et al. (i.e., 
solving system (5)) would not produce meaningful results. 
However, as we show in this section, the matrix C still ex¬ 
hibits a particular structure which can be exploited to drive 
the matching process. 

We assume to be given a full shape M and a part thereof 
J\T C AT. We further denote by J\T = M \N the remaining 
vertices of AT. The manifolds M and N are discretized as 
triangular meshes with m and n vertices, respectively, and 
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h = m — n. The scenario we consider in this paper con¬ 
cerns the problem of matching an approximately isometric 
deformation of part A f to the full shape M (part-to-whole 
matching). Our goal is to characterize the eigenvalues and 
eigenvectors of the Laplacian in terms of perturbations 
of the eigenvalues and eigenvectors of the Laplacians L jq- 
and [MH88]. We tacitly assume that homogeneous Neu¬ 
mann boundary conditions (2) apply. 


3.1. Block-diagonal case 

For the simplicity of analysis, let us first consider a sim¬ 
plified scenario in which Af and Af are disconnected , i.e., 
there exist no links between the respective boundaries dAf 
and dAf. W.l.o.g., we can assume that the vertices in A4 
are ordered such that the vertices in Af come before those 
in Af. With this ordering, the m x m Laplacian matrix Lj^ 
is block-diagonal, with an n x n block Ljq and an h x n 
block L A t- The (sorted) eigenvalues of form a mixed 
sequence composed of the eigenvalues from Ljq and Ljq. 
Similarly, the eigenvectors of L m correspond to the eigen¬ 
vectors of the two sub-matrices, zero-padded to the correct 
size (Fig. 3). 

Structure of C under partiality. Suppose we are given the 
first k Laplace-Beltrami eigenvalues of the full shape A4 and 
of its part Af. Since the spectrum of is an interleaved 
sequence of the eigenvalues of L jq and L^, only the first 
r < k eigenvalues of Ljq will appear among the first k eigen¬ 
values of L. The remaining k — r eigenvalues of Ljv will 
only appear further along the spectrum of L(see Fig. 6 
for an example where k = 50 and r = 21). The same argu¬ 
ment holds for the associated eigenfunctions, as illustrated 
in Fig. 4: if <|>; is an eigenvector of Ljq, then also has an 
eigenvector \|fj such that = T\|/y, where T = (l n xn, 0) T 
and i < j. 

This analysis leads us to the following simple observation: 
the partial functional map between Af and A4 is represented 



Figure 3: The eigenvalues and eigenvectors of a block- 
diagonal Laplacian L j^i are an interleaved sequence of the 
eigenpairs from the two blocks L jq and ~Ljq. 


in the spectral domain by the matrix of inner products Cij = 
(T<| >i,Vj)M, which has a slanted-diagonal structure with a 
slope r/k (see examples in Figs. 1, 4 where this structure 
is manifested approximately). Consequently, the last k — r 
columns of matrix C are zero such that r = rank(C). The 
value r can be estimated by simply comparing the spectra of 
the two shapes, as shown in Fig. 6. Note that this behavior is 
consistent with Weyl’s asymptotic law [Weyll], according 
to which the Laplacian eigenvalues grow linearly, with rate 
inversely proportional to surface area. 


3.2. Perturbation analysis 

We will now show that these properties still approximately 
hold when the Laplacian matrix Lis not perfectly block- 
diagonal, i.e., when Af and Af are joined along their bound¬ 
aries. Roughly speaking, the main observation is that in this 
case as well the matrix C has a slanted diagonal structure, 
where the diagonal angle depends on the relative area of the 
part, and the diagonal ‘sharpness’ depends on the position 
and length of the cut. 

Here we assume w.l.o.g. that within Af the boundary ver¬ 
tices dAf are indexed at the end, while within Af the bound¬ 
ary vertices dAf are indexed at the beginning. Then, there is 
a boundary band B = dAf U dAf such that only the entries 
of the Laplacians L jq and between vertices in B are af¬ 

fected by the cut (Fig. 5). 


We define the parametric matrix 

L (0 = 


L M 

0 N 

\ + ,( P -V 

p 

0 

L A T , 

) + V T 

P AT 


where 

Pw = 


0 0 

0 Djv" 


la r 


= (^ 0 


0 0 


, p = 


(7) 


0 0 
E 0 


are matrices of size nxn,hxh, and nxh, respectively. Here 
D/v" and represent the variations of the Laplacians L/v 
and L A T within nodes in dAf and dAf respectively, while E 
represents the variations across the boundary. The parame¬ 
ter t is such that L(l) = while for t = 0 we get back 



Figure 5: The matrix L (t) is obtained as a perturbation of 
the block-diagonal Laplacian in the boundary band (shown 
in green). 


submitted to COMPUTER GRAPHICS Forum (12/2015). 










































E. Rodola et al. / Partial Functional Correspondence 


5 




4>2 fa fa fa < i >6 < i >7 <|>8 fa 



Figure 4: First ten eigenfunctions of a full shape M and two parts A/i,A /*2 C M with different surface area. All eigenfunc¬ 
tions of the partial shapes have a corresponding eigenfunction \|/; on the full shape for some i; the correspondence between 
eigenfunctions follows from the correspondence between eigenvalues (see also Fig. 6). This is reflected in functional maps with 
different diagonal slopes, where the slope depends on the area ratios of the two surfaces (by WeyVs law). 


to the disconnected case of Fig. 3. In what follows, we per¬ 
form a differential analysis at t = 0 and therefore analyze the 
change in eigenvalues and eigenvectors of L(?) as we inter¬ 
polate between the seen (A f) and unseen (J\f) parts for t = 0 , 
to the full shape M for t = 1 . 

Note that with the appropriate ordering of the vertices, the 
matrices D/v and will have a band-diagonal structure. 
In fact, a cut through an edge will affect the values of the 
discrete Laplacian matrix L only at the entries corresponding 
to the vertices at the extremities of the edge, and to the edges 
laying in the same triangle as the cut edge. For example, 
looking at Fig. 2, a cut through edge (/, j) will affect the 
diagonal entries la and Ijj as well as the off-diagonal entries 
lih , ljh, and Ijk. Note also that the continuity of the cut 
implies that two of the four off-diagonal entries will be cut 
as well, leaving no more than two affected edges on any side 
of the cut. As a result, the entries of the Laplacian affected 
by the cut correspond to the nodes and edges in a path along 
the boundary of the cut. 

We further note that although here we consider the cotan¬ 
gent Laplacian for simplicity of analysis, similar results hold 
for Laplacians that are not strictly local, but locally dom¬ 
inant [BSW09b], i.e. y most of their L 2 norm is due to the 
elements in a tight boundary layer. 

Theorem 1 LetL/v + rP^ = &(t) T A(t)&(t), where A (t) = 
diag(A,i (t),...,X n (t)) is a diagonal matrix of eigenvalues, 
and 0(t) are the corresponding eigenvectors. The derivative 
of the non-trivial eigenvalues is given by 

“77 A,/ — ^ (J*J\f)vwtyivtyiw = P/V^r ( 8 ) 

v,wEdAf 


Theorem 1 establishes that the (first-order) change in the 
eigenvalues of the partial shape A f only depends on the 
change in the Dirichlet energy of the corresponding eigen¬ 
vectors along the boundary dj\f (recall that for eigenvector 
<|> i , the Dirichlet energy is defined as (L yf + jPa/')^- = 

Pa/"^/)- This means that the eigenvalues are 
perturbed depending on the length and position of the cut. 
Note that the estimate given in Eq. ( 8 ) can not typically be 
computed directly, as this would assume knowledge of the 
correspondence between J\f and AT. However, by virtue of 
this result, we can establish approximate correspondence be¬ 
tween the eigenvalues Xj^ of L/y/ and a subset of the eigen¬ 
values X^ of L^/i (which are now not exactly equal as in 
the block-diagonal case). We do this in order to estimate the 
slope of C. Specifically, we compute 

r = max{/1 Xf^ < max }; (9) 

the slope of C can now be estimated as r/k, as explained in 
Section 3.1 (see also Fig. 6 for a visual illustration of the 
estimation of r). 


Theorem 2 Assume that Ljv has distinct eigenvalues (A; 7 ^ 
Xj for i f j ), and furthermore, the non-zero eigenvalues 
are all distinct from the eigenvalues of (A* 7 ^ Xj for 
all /, j). Let La/' + t~Pj\f = &(t) T A(t)&(t), where A (t) = 
diag(A,i(f),..., X n (t)) is a diagonal matrix of eigenvalues, 
and &(t) are the corresponding eigenvectors. Then, the 
derivative of the non-constant eigenvector is given by 


d n 

Tt+i = E 


dt 


j= 1 
i¥i 




(10) 


Proof: See Appendix C. 


Proof: See Appendix C. 
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Figure 6: Neumann spectra of a full shape and a part of 
it. The eigenvalues of the partial shape (in red) are approx¬ 
imately preserved under the partiality transformation (see 
Theorem 1 ), and appear perturbed in the spectrum of the 
full shape (in blue). This simple observation allows us to es¬ 
timate the diagonal slope of the functional map relating the 
two shapes; in this example, the slope is equal to 21/50. 



Remark. If shares some eigenvalues with L/y/, the sec¬ 

ond sum in (10) would be slightly different [MH 88 ], but 
would still only have support over M. 

We conclude from Theorem 2 that the perturbation associ¬ 
ated with the partiality transformation gives rise to a mixing 
of eigenspaces. The second summation in (10) has support 
over M and thus provides the completion of the eigenfunc¬ 
tion on the missing part. The first summation in (10) is re¬ 
sponsible for the modifications of the eigenvectors over the 
nodes in J\f. Here the numerator has a term tyj P 'N§j which, 
since Djq- is band-diagonal and diagonally dominant, acts as 
a dot product of the eigenvectors over the boundary band. 
This points to large mixing of eigenvectors with a strong co¬ 
presence near the boundary. In turn, the term X[ — Xj at the 
denominator forces a strong mixing of eigenvectors corre¬ 
sponding to similar eigenvalues. This results in an amplifi¬ 
cation of the variation for higher eigenvalues, as eigenvalues 
tend to densify on the higher end of the spectrum, and ex¬ 
plains the funnel-shaped spread of the matrix C visible at 
high frequencies (see Fig. 8 ). 

Similarly to the case of eigenvalues, the eigenvectors are 
also perturbed depending on the length and position of the 
cut. The variation of the eigenvectors due to the mixing 
within the partial shape can be reduced either by shortening 
the boundary of the cut, or by reducing the strength of the 
boundary interaction. The latter can be achieved by select¬ 
ing a boundary along which eigenvectors with similar eigen¬ 
values are either orthogonal, or both small. The boundary 
interaction strength can be quantified by considering the fol¬ 
lowing function (we refer to Appendix C for a derivation): 



Figure 7: Left: A model is cut in two different ways (red 
and green curves) with cuts of same length. The off-diagonal 
dispersion depends mainly on the position of each cut. Func¬ 
tion f (11) is plotted over the model. Middle: Ground-truth 
functional map between the complete model and the partial 
shape produced by the red cut (top), and values off along 
the cut (bottom). Right: Plots associated to the green cut. 


/(v) = 


I 

i,j= 1 
H 



2 


( 11 ) 


Fig. 7 shows an example of two different cuts with different 
interaction strengths, where the function / is plotted on top 
of the cat model. The cuts plotted in the figure have the same 
length, but one cut goes along a symmetry axis of the shape 
and through low values of /, while the other goes through 
rather high values of /. This is manifested in the dispersion 
of the slanted diagonal structure of the matrix C (larger in 
the second case). 


4. Partial functional maps 

As stated before, throughout the paper we consider the set¬ 
ting where we are given a full model shape M and another 
query shape J\f that corresponds to an approximately isomet- 
rically deformed part M r C A4. 

Following [BB08], we model the part M 1 by means of 
an indicator function v : M {0, 1 } such that v(x) = 1 if 
x C M r and zero otherwise. Assuming that v is known, the 
partial functional correspondence between A f and M can be 
expressed as T f — vg, where v can be regarded as a kind of 
mask, and anything outside the region where v = 1 should 
be ignored. Expressed w.r.t. bases {(()/}/>i and {\|//}i>l> 
the partial functional correspondence takes the form CA = 
B(v), where B(v) denotes a matrix of weighted inner prod¬ 
ucts with elements given by bij(v) = J M v(x)\]fi(x)gj(x)dx 
(when v(x) = 1, B is simply the matrix of Fourier coeffi¬ 
cients defined in (4)). 

This brings us to the problem we are considering through¬ 
out this paper, involving optimization w.r.t. correspondence 
(encoded by the coefficients C) and the part v, 

min ||CA-B(r|(v))|| 2 ,l +Pcorr(C) + Ppart(v), (12) 

C,v 

where r\(t) = \ (tanh( 2 r — 1 ) + 1 ) saturates the part indica- 
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tor function between zero and one (see below). Here p CO rr 
and Ppart denote regularization terms for the correspondence 
and the part, respectively; these terms are explained below. 
We use the matrix norm (equal to the sum of L 2 -norms 
of matrix columns) to handle possible outliers in the corre¬ 
sponding data, as such a norm promotes column-sparse ma¬ 
trices. A similar norm was adopted in [HWG14] to handle 
spurious maps in shape collections. 

Note that in order to avoid a combinatorial optimization 
over binary-valued v, we use a continuous v with values in 
the range (—00, +00), saturated by the non-linearity rj. This 
way, r| (v) becomes a soft membership function with values 
in the range [ 0 , 1 ]. 

Part regularization. Similarly to [BB08,PBB 13], we try to 
find the part with area closest to that of the query and with 
shortest boundary. This can be expressed as 

Ppart (v) = p\ ^area(A0 - J^x\{v)dx^ (13) 

+ P 2 [ $(v)\\V M v\\dx, 

JM 

where £(f) ^ 3 (r\{t) — and the norm is on the tangent 
space. The /^-term in (13) is an intrinsic version of the 
Mumford-Shah functional [MS 89], measuring the length of 
the boundary of a part represented by a (soft) membership 
function. This functional was used previously in image seg¬ 
mentation applications [VC02]. 

Correspondence regularization. For the correspondence, 
we use the penalty 

Pcorr(C) = W ||CoW||g+/4£(C T C$ 

ifj 

+ W £((C T C ) u -di) 2 , (14) 

i 

where o denotes Hadamard (element-wise) matrix product. 
The ^ 3 -term models the special slanted-diagonal structure of 
C that we observe in partial matching problems (see Fig. 8 ); 
the theoretical motivation for this behavior was presented 
in Sec. 3. Here, W is a weight matrix with zeros along the 
slanted diagonal and large values outside (see Fig. 8 ; details 
on the computation of W are provided in Appendix A). 

The /^ 4 -term promotes orthogonality of C by penalizing 
the off-diagonal elements of C T C. The reason is that for 
isometric shapes, the functional map is volume-preserving, 
and this is manifested in orthogonal C [OBCS* 12]. Note that 
differently from the classical case (i.e., full shapes), in our 
setting we can only require area preservation going in the 
direction from partial to complete model, as also expressed 
by the p i-term in (13). For this reason, we do not impose 
any restrictions on CC T and we say that the matrix is semi- 
orthogonal. 

Finally, note that due to the low-rank nature of C we can 



Figure 8: A partial functional map Cfrom N to M has a 
slanted-diagonal structure (second row, left). The low-rank 
nature of such a map is also manifested in its singular val¬ 
ues (bottom right). If the map is volume-preserving, then its 
full-rank sub-matrix is orthogonal: observe how C T C ap¬ 
proximates the identity, with a trail of small values along the 
diagonal corresponding to the almost-zero block of C. 


not expect the product C T C to be full rank. Indeed, we ex¬ 
pect elements off the slanted diagonal of C to be close to 

zero and thus C T C « ^ J . The ^ 5 -term in (14) models 

this behavior, where vector d = (d\^... : df) determines how 
many singular values of C are non-zero (the estimation of d 
is straightforward, and described in Appendix A). 

Remark. The fact that matrix C is low-rank is a direct con¬ 
sequence of partiality. This can be understood by recalling 
from Eq. (4) that the (non-truncated) functional map repre¬ 
sentation amounts to an orthogonal change of basis; since in 
the standard basis the correspondence matrix is low-rank (as 
it contains zero-sum rows), this property is preserved by the 
change of basis. 

In Fig. 8 we show an example of a ground-truth partial func¬ 
tional map C, illustrating its main properties. 

Alternating scheme. To solve the optimization prob¬ 
lem (12), we perform an alternating optimization w.r.t. to C 
and v, repeating the following steps until convergence: 

C-step: Fix v*, solve for correspondence C 

nun ||CA-B(t|(v*))|| 2 ,i +Pcorr(C). (15) 
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Figure 9: Aw example of the matching process operating on two shapes from TOSCA. The algorithm alternatingly optimizes 
over corresponding part (top row) and functional correspondence (bottom row). Corresponding points between full and partial 
shape are shown with the same color. This solution was obtained by using 30 eigenfunctions on both manifolds. 


V-step: Fix C*, solve for part v 

min ||C*A-B(r|(v))||2 ) i +p P art(v). (16) 

A practical example of the alternating scheme applied to a 
pair of shapes is shown in Fig. 9. 


5. Implementation 

We refer to Appendix A for the discretization of the regular¬ 
ization terms appearing in (13) and (14). 

Numerical optimization. We implemented our matching 
framework in Matlab/C++ using the manifold optimization 
toolbox [BMAS14]. Each optimization step was performed 
by the method of nonlinear conjugate gradients. Detailed 
derivations of the involved gradients can be found in Ap¬ 
pendix B. We initialize the alternating scheme by fixing 
v* = 1 (a vector of m ones), C = W, and by optimizing over 
C. In all our experiments we observed convergence in 3-5 
outer iterations (around 5 mins, for a pair of shapes). 

Refinement. In order to account for noisy data, we also run 
a refinement step after each C-step. Specifically, assume C* 
is a local optimum of problem (15), and consider the term 
||C*<I> T -'P T n|| f( where n is a left-stochastic binary ma¬ 
trix assigning each column of to the nearest column of 
C*<I> T ; this is done by n nearest-neighbor searches in R k , 
one per column. Given the optimal II*, we solve for the 
map C minimizing ||C<I> T — 'P T II* \\p plus the / 14 , p$ terms 
of Eq. (14). We alternate the C* and II* steps until conver¬ 
gence. This refinement step can be seen as a generalization to 
partial maps of the ICP-like technique found in [OBCS* 12], 
and can be interpreted as an attempt to improve the align¬ 
ment between the spectral embeddings of the two shapes. 
Further note that matrix II* encodes the point-wise corre¬ 
spondence between J\f and A4, which is used to evaluate the 
accuracy of our method. 


6. Experimental results 

Datasets. As base models, we use shapes from the TOSCA 
dataset [BBK08], consisting of 76 nearly-isometric shapes 
subdivided into 8 classes. Each class comes with a “null” 
shape in a standard pose (extrinsically bilaterally symmet¬ 
ric), and ground-truth correspondences are provided for all 
shapes within the same class. In order to make the datasets 
more challenging and avoid compatible triangulations, all 
shapes were remeshed to 10K vertices by iterative pair con¬ 
tractions [GH97]. Then, missing parts were introduced in the 
following ways:^ 

Regular cuts. The null shape of each class was cut with a 
plane at 6 different orientations, including an exact cut along 
the symmetry plane. The six cuts were then transferred to 
the remaining poses using the ground-truth correspondence, 
resulting in 456 partial shapes in total. Some examples are 
shown in Fig. 6 and 8 . 

Irregular holes. Given a shape and an “area budget” de¬ 
termining the fraction of area to keep (40%, 70%, and 90%), 
we produced additional shapes by an erosion process applied 
to the surface. Specifically, seed holes were placed at 5, 25, 
and 50 farthest samples over the shape; the holes were then 
enlarged to meet the specified area budget. The total num¬ 
ber of shapes produced this way was 684. Examples of this 
dataset are shown in Fig. 11 and 13. 

Range images. We simulated range images by taking or¬ 
thographic projections of the original TOSCA shapes from 
different viewpoints. Each range image was produced via ray 
casting from a regular grid with a resolution of 100 x 150 
pixels. Examples are shown in Fig. 13. 

Point clouds. Point clouds were generated by taking a sub¬ 
set of shapes from the first two datasets. Each partial shape 


t The datasets together with code for our method are avail¬ 
able for download at http://vision.in.tum.de/data/ 
datasets/partial. 
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Figure 10: Correspondence quality of different methods 
evaluated using the Princeton protocol on partial TOSCA 
shapes with regular cuts (solid) and irregular holes (dotted). 


was then resampled uniformly to 1000 farthest points, and 
the tessellation removed. See Fig. 13 for examples. 

Where not specified otherwise, we use 120 random partial 
shapes for the first dataset and 80 for the second, equally 
distributed among the different classes. Each partial shape is 
then matched to the null shape of the corresponding class. 


Error measure. For the evaluation of the correspondence 
quality, we used the Princeton benchmark protocol [KLF1 1] 
for point-wise maps. Assume that a correspondence algo¬ 
rithm produces a pair of points (x,y) eNfxM, whereas the 
ground-truth correspondence is (x,y*). Then, the inaccuracy 
of the correspondence is measured as 


e(x) 


djvdyy) 

area(.M) 1 / 2 ’ 


(17) 


and has units of normalized length on M (ideally, zero). 
Here dy^ is the geodesic distance on M. The value e(jc) 
is averaged over all shapes J\f. We plot cumulative curves 
showing the percent of matches which have error smaller 
than a variable threshold. 


Methods. We compared the proposed method with (full) 
functional maps [OBCS*12], elastic net [RTH*13], and the 
voting method of [SY14] using the code provided by the re¬ 
spective authors. 

Local descriptors. Due to the particular nature of the prob¬ 
lem, in all our experiments we only make use of dense, local 
descriptors as a data term. This is in contrast with the more 
common scenario in which full shapes are being matched 
- thus allowing to employ more robust, globally-aware fea¬ 
tures such as landmark matches, repeatable surface regions, 
and various spectral quantities [OBCS*12]. In our experi¬ 
ments, we used the extrinsic SHOT [TSDS10] descriptor, 
computed using 10 normal bins (352 dimensions in total). 




Figure 11: Correspondence quality (in terms of mean 
geodesic error, in % of diameter) obtained by different 
methods at increasing levels of partiality. Other methods 
show significant performance drop with increasing partial¬ 
ity, while the performance of our method is nearly-constant. 


As opposed to [ART15, PBB13] which ignore points close 
to the boundary in order to avoid boundary effects, in our 
formulation we retain all shape points. 


6.1. Sensitivity analysis 

We conducted a set of experiments aimed at evaluating the 
sensitivity of our approach to different parametrizations. In 
order to reduce overfitting we only used a subset of TOSCA 
(regular cuts), namely composed of the cat and victoria 
shape classes (20 pairs). 



Figure 12: Correspondence quality obtained on a subset of 
TOSCA at increasing rank (reported as labels on top of the 
curves). Note the opposite behavior of the baseline approach 
and our regularized partial matching. 
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Figure 13: Examples of partial functional correspondence obtained with our method on meshes and point clouds from the 
proposed datasets. Notice how regions close to the boundary are still accurately matched despite the noisy descriptors. 


Rank. In the first experiment we study the change in ac¬ 
curacy as the rank of the functional map is increased; this 
corresponds to using an increasing number of basis func¬ 
tions for the two shapes being matched. For this experi¬ 
ment we compare with the baseline method of Ovsjanikov et 
al. [OBCS* 12] by using the same dense descriptors as ours. 
For fair comparisons, we did not impose map orthogonal¬ 
ity or operator commutativity constraints [OBCS* 12], which 
cannot obviously be satisfied due to partiality. The results 
of this experiment are reported in Fig. 12. As we can see 
from the plots, our method allows to obtain more accurate 
solutions as the rank increases, while an opposite behavior 
is observed for the other method. 

Representation. Our method is general enough to be ap¬ 
plied to different shape representations, as long as a proper 
discretization of the Laplace operator is available. In Fig. 13 
we show some qualitative examples of correspondences pro¬ 
duced by our algorithm on simulated point clouds and depth 
maps. Here we use the method described in [BSW09a] to 


construct a discrete Laplacian on the point clouds. This is 
traditionally considered a particularly challenging problem 
in robotics and vision applications, with few methods cur¬ 
rently capable of giving satisfactory solutions without ex¬ 
ploiting controlled conditions or domain-specific informa¬ 
tion ( e.g ., the knowledge that the shape being matched is that 
of a human). These are, to the best of our knowledge, the best 
results to be published so far for this class of problems. 


6.2. Comparisons 

We compared our method on the cuts and holes datasets 
(200 shape pairs in total) ; the results are shown in Fig. 10. 
As an additional experiment, we ran comparisons against 
[OBCS* 12] across increasing amounts of partiality. The ra¬ 
tionale behind this experiment is to show that, at little or 
no partiality, our approach converges to the one described 
in [OBCS* 12], currently among the state of the art in non- 
rigid shape matching. However, as partiality increases so 
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does the sensitivity of the latter method. Fig. 11 shows the 
results of this experiment. 

Parameters for our method were chosen on the basis of 
the sensitivity analysis. Specifically, we used k = 100 eigen¬ 
functions per shape, and set p\ = = 1 , ^4 = p$ = 10 3 , and 

P2 = 10 2 - The different orders of magnitude for the p coeffi¬ 
cients are due to the fact that the regularizing terms operate 
at different scales. We also experimented with other values, 
but in all our tests we did not observe significant changes in 
accuracy. Additional examples of partial matchings obtained 
with our method are shown in Fig. 13. 

7. Discussion and conclusions 

In this paper we tackled the problem of dense matching of 
deformable shapes under partiality transformations. We cast 
our formulation within the framework of functional maps, 
which we adapted and extended to deal with this more chal¬ 
lenging scenario. Our approach is fully automatic and makes 
exclusive use of dense local features as a measure of simi¬ 
larity. Coupled with a robust prior on the functional corre¬ 
spondence derived from a perturbation analysis of the shape 
Laplacians, this allowed us to devise an effective optimiza¬ 
tion process with remarkable results on very challenging 
cases. In addition to our framework for partial functional 
correspondence we also introduced two new datasets com¬ 
prising hundreds of shapes, which we hope will foster fur¬ 
ther research on this challenging problem. 

One of the main issues of our method concerns the exis¬ 
tence of multiple optima, which is in turn related to the pres¬ 
ence of non-trivial self-isometries on the considered man¬ 
ifolds. Since most natural shapes are endowed with intrin¬ 
sic symmetries, one may leverage this knowledge in or¬ 
der to avoid inconsistent matchings. For example, includ¬ 
ing a smoothness prior on the correspondence might allevi¬ 
ate such imperfections and thus provide better-behaved solu¬ 
tions. Secondly, since the main focus of this paper is on tack¬ 
ling partiality rather than general deformations, our current 
formulation does not explicitly address the cases of topo¬ 
logical changes and inter-class similarity (e.g., matching a 
man to a gorilla). However, the method can be easily ex¬ 
tended to employ more robust descriptors such as pairwise 
features [RABT13, vKZH13], or to simultaneously optimize 
over ad-hoc functional bases on top of the correspondence. 
Finally, extending our approach to tackle entire shape col¬ 
lections, as opposed to individual pairs of shapes, represents 
a further exciting direction of research. 
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the sake of compactness). Detailed gradients for each term are given 
in Appendix B. 


Mesh parametrization. We denote by S a triangle mesh of 
n points, composed of triangles Sj for j = 1,.. . ,m. In the fol¬ 
lowing derivations, we will consider the classical triangle-based 
parametrization described by the charts xj : M 2 —M 3 

(a, P) = xja + a{xj , 2 ~Xj, i) + P(y /, 3 -xj, i), (18) 

with a E [0,1] and p E [0,1 — a]. With Xj^ E M 3 we denote the 3D 
coordinates of vertex k E {1,2,3} in triangle Sj. 

Each triangle Sj is equipped with a discrete metric tensor with 
coefficients 



where Ej = \\xj 2 — x/.i || 2 , Fj — (xjp — Xjy,Xj^ — Xj _\), and Gj = 
| |xy _3 — Xj_ 1 11 2 . The volume element for the y-th triangle is then given 
by = y/EjGj-Fj. 


Integral of a scalar function. Scalar functions / : S M are 

assumed to behave linearly within each triangle. Hence, f(x( a, p)) 
is a linear function of (a, P) and it is uniquely determined by its 
values at the vertices of the triangle. The integral of / over Sj is 
then simply given by: 


^‘^ 1 ““/(a,p) v / det^^a (20) 

= j f f(0,0)(1 — a — P) +/(l,0)a + /(0,1)P y/detgydpda 
= 1 (/(0,0) + /(1,0) + f(0, 1)) 

( 21 ) 


= 3 (f(0,0) + /(1,0) +/(0, l))area(S;) 


where/( 0 , 0 ) =/(x y ,i),/(l, 0 ) =/(x ; - 2 ),and/( 0 , 1 ) =/(x ; - 3 ). 


Gradient of a scalar function. For the intrinsic gradient of / 
we get the classical expression in local coordinates: 



where we write f a to denote the partial derivative = fjp — f j. \ 
and similarly for /p. The norm of the intrinsic gradient over triangle 
Sj is then given by: 


l|V/|| (23) 

= VWINT) 



h ) 



fa 

h 



h ) 



fa A 1 

h ) \/det Yj 




det gj 


(24) 


Appendix A - Discretization 

We describe more in depth the discretization steps and give the im¬ 
plementation details of our algorithm (we skip trivial derivations for 
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Note that, since we take / to be linear, the gradient V/ is constant 
within each triangle. We can then integrate V/ over Sj as follows: 


[ \\Vf(x)\\^MgjdadV 

JSi 


(25) 


L 


: j\ 


fjfr 2,/ a ,/ 3 /) • f>Ej 

det gj 


^/delgjdaxlfi 


= [ Jf&Gj-lfafiFj + flEj dad$ 

JSj 

= \y/f&G J -2f a Wi+tfEj. 


(26) 


In the following, we write JV and A4 to denote the partial and 
full shape respectively. Further, let be the first k eigen¬ 
values of the Laplacian on M, and similarly for The 

functional map C has size k x k. 


Mumford-Shah functional (// 2 -term). Following Equations 
( 21 ) and (26), we immediately obtain: 

jk(v)\\Vv\\dx 

m p 

= E / Ua,mVv\\y^ijdpda 

j=\ JS J 

m t - r 

= L V v a G J - 2 v a v V F j + vjEj / B,(a, fydfida 

i m - 

« g £ y/vkGj - 2v„vpFy + ^£y(5(0,0) + §(1,0) + 5(0,1)), 

where £( 0 , 0 ) = £(v(x ; -,i))> £( 1 , 0 ) = £(v(x/, 2 )), and £( 0 , 1 ) = 

U v ( x j, 3))- 

Weight matrix ( 1/3 -term). Recall from Section 3 and Figure 6 
that an estimate for the rank of C can be easily computed as 

r = max{/ | < max^-f 4 } . (27) 

j 

We use this information in order to construct the weight matrix W, 
whose diagonal slope directly depends on r. 

To this end, we model W as a regular k x k grid in M 2 . The slanted 
diagonal of W is a line segment 5(t) = p +1 j^y with t E M, where 
p = ( 1 , 1 ) T is the matrix origin, and n = (l,r/k) T is the line di¬ 
rection with slope r/k. The high-frequency spread in C is further 
accounted for by funnel-shaping W along the slanted diagonal. We 
arrive at the following expression for W: 

Wij = JL x ((i,j) T -p)|| , (28) 

where the second factor is the distance from the slanted diagonal 5, 
and a E M+ regulates the spread around 5. In our experiments we 
set a = 0.03. 


Orthogonality (^ 4,^5 -terms). For practical reasons, we incor¬ 
porate the off-diagonal and diagonal terms within one term with 
the single coefficient ^ 5 . In addition, we rewrite the off-diagonal 
penalty using the following equivalent expression: 

L(C T C)?. = ||C T C||p-£(C T C)|. (29) 

Hi i 


Vector d E is constructed by setting the first r elements (accord¬ 
ing to (27)) equal to 1, and the remaining k — r elements equal to 
0 . 


Appendix B - Gradients 

We find local solutions to each optimization problem by the (nonlin¬ 
ear) conjugate gradient method. In this Section we give the detailed 
gradient derivations of all terms involved in the optimization. 

In order to keep the derivations practical, we will model function 
v by its corresponding n-dimensional vector v. Note that, depending 
on the optimization step, the gradients are computed with respect to 
either v or C. 


Data term (w.r.t. v). Let q denote the number of correspond¬ 
ing functions between the two shapes. Matrices F and G contain 
the column-stacked functions defined over M and Ai and have 
size nx q. The respective projections onto the corresponding func¬ 
tional spaces 4> and are stored in the k x q matrices A and B 
respectively. Let us write Hj 7 to identify the elements of the matrix 
CA — B(r|(v)), we then have 


— ||CA-B(n(v))|| 2 ,i 

ov„ 




■E Mi 1 " 11 


j= 1 v=i / i= 1 9v p 

r'n iTj I Y 


Since B(r|(v))y = we have: 


dv„ Kii - dv„ [CAiJ ByI “ ^ Fw ^ T,(v ' ) ■ 


Finally: 


-—r| (v p ) - 1 - tanh 2 ( 2 v p - 1 ). 

dVn 


(30) 


(31) 


(32) 


Area term (//i-term w.r.t. v). The derivative of the discretized 
area term is: 

— 2 ^ V Av) 1 _ y Am Sm n( v ’/’) 

where (Sm); and (Sjq-)i are the local area elements associated with 
the z'-th vertex of meshes Ai and AT respectively. For the derivative 
of t|(vp) see equation (32). 


Mumford-Shah functional (t/ 2 -term w.r.t. v). Computing 
the gradient V v / 5 ^(v) II Vv||dx involves computing partial deriva¬ 
tives of ^(v) with respect to v. These are simply given by: 



tanh(2v£ — 1) 
4 o 2 


1 — tanh 2 ( 2 vfc — 1 ) - 
2 a 2 ^ 


tanh(2v£ — 1) 

4 o 2 


submitted to COMPUTER GRAPHICS Forum (12/2015). 















14 


E. Rodola et al. / Partial Functional Correspondence 


In the following derivations we set Dj = ^v^Gj — 2v a v$Fj + v^Ej, 
and Dj = 0 whenever Vv = 0. The gradient of the Mumford-Shah 
functional is then composed of the partial derivatives: 

[ ^»l|Vv||<fc 
ov k Js 

m ? r 

= ^ v )H Vv ll 

pi °n Js, 

= 1 I ^-Dj&M + Uvj, 2 ) + 5 (vy, 3 )) 

6 jeN(k ) 

= J I ffiv i )+?to)+^ 3 ))/-D / +D i ^(n) 

6 ^ 9v* 

= ' E (5W+5(v^) + «v W ))rf-^-Dj+D y3 l5(v t ) 

6 ye V« 2Djdv k dv k 

= ; I (5(v.) + 5(vi )2 )+fe))^ + Oi^(v t ) 

6 j£N(k) dv k 

where we write ^ = ^((v* - v j,2)( G j - Fj) + (v k - v jj3 )(Ej - 
Fj)), and j G A(k) are the indices of the triangles containing the 
k -th vertex. Note that we slightly abuse notation by writing v k , Vj t 2 , 
and v J; 3 to denote the three vertices of the y'-th triangle, even though 
in general the ordering might be different depending on the triangle. 

Data term (w.r.t. C). The derivative of the data term with respect 
to C is similar to (30). The only difference is in the partial derivative: 


d d * 

,)c H " “ ;ir E (: «V; - 


9Cn 


C pq A qj if i = p 


0 


otherwise. 


(33) 


Weight matrix (p 3 -term w.r.t. C). This is simply given by: 

ac 


^-||CoW|| 2 =C p? (W p? ) 2 . 


(34) 


and A(f ) is diagonal and &(t) orthogonal. 

Following [MH 88 ], if all the eigenvalues are distinct, then we can 
compute the derivatives of the eigenvalues at t = 0 as 

K = <t>y A> ( (37) 

where A', the derivative of A (t), and the eigenvectors <J> ; are consid¬ 
ered at t = 0. In fact, differentiating (35), we obtain 

A 4> + A<S>' = <t>'A + 0>A' . (38) 

Left-multiplying both sides by <I> T , setting <t> ; = 4>B for a matrix B 
to be determined, and recalling that 4> T A4> = A, we have 

4> t A 4> + AB = BA + A' , (39) 


from which 

diag(A / ) = diag(<I> T A / 4>) + diag(AB - BA) = diag(4> T A'0>). 

(40) 

Going back to our case, we take the simplifying assumptions that 
Lyvi and L/v" do not have repeated eigenvalues. Let 

L/v" = 4 > t A4> (41) 

L m = <PA<P (42) 

be the spectral decompositions of Ly^ and respectively. Accord¬ 
ing to the previous result, we can write the derivative of A; eigen¬ 
value ofLyy and, thus, of L(0), as: 

X| = fTp^= £ (Pa/-M;v<Iw- (43) 

v,w£d Af 

Theorem 2 Assume that Ljy has distinct eigenvalues (A* 7 ^ Ay for 
i 7 ^ j), and furthermore, the non-zero eigenvalues are all distinct 
from the eigenvalues of L^ (A/ 7 ^ Ay for all i, j). Let L yy + tPj\f = 
4>(r) T A(r)4>(r), where A (t) = diag(Ai(f),. .., X n (t)) is a diagonal 
matrix of eigenvalues, and 4>(r) are the corresponding eigenvectors. 
Then, the derivative of the non-constant eigenvector is given by 


Orthogonality (jli 4 5 -term w.r.t. C). The gradient of the last 
term can be finally obtained as: 

C T C§- E(C T C)i + £((C T C)„ - dif 



j¥i 


(44) 


Proof: Under the same assumptions as for the previous theorem, 
from (39) we have 


= 4(CC r C) p? + 2£ 

i 

= 4[(CC r C)p 9 — dqCpql . 


-EC li + (Eel,. - d t ) -y Ec* 2 , 

k a ^pq k k °^Pq k 


Appendix C - Perturbation Analysis 

Theorem 1 Let Ly^ + rPyy — <I>(f) T A(r)<J>(r), where A(f) = 
diag(Ai(f),..., X n (t)) is a diagonal matrix of eigenvalues, and &(t) 
are the corresponding eigenvectors. The derivative of the non-trivial 
eigenvalues is given by 

(35) 

dt v,w£dAT 

Proof: Let A(f) be a symmetric real n x n matrix parametrized by 
t CT CM, with &(t) and A (t) being the eigenvector and eigenvalue 
matrices, i.e ., for all t G T we have 

A(0*(t) = *MA(0 (36) 


(* T A'*) w + (AB)y = (BA)y + (A')y 

♦fA '^j + Xibij = bijXj+0, 

from which 

lJ Xj-Xi ’ 


(45) 

(46) 


(47) 


due to the orthogonality of 4>(t) we have that B is skew-symmetric, 
and thus, bn = 0. From the relation O' = 4>B we obtain 


d 

dt 


< t>, = Em,- = E 

j¥* j¥* 



(48) 


Going back to our case, recall that the set of eigenvalues of L(0) 
is the union of the eigenvalues of Lyy and and the correspond¬ 
ing eigenvectors are obtained from those of Lyy and Lyy by padding 
with zeros on the missing parts. Denoting with A; and <!>, the eigen¬ 
values and corresponding eigenvectors of L yy, and Xj and <{) / the 
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eigenvalues and corresponding eigenvectors of L^, we have 

d " *[P f <I>, T P^ X 

j¥i 


(49) 


Boundary interaction strength. We can measure the variation 
of the eigenbasis as a function of the boundary B splitting M. into 
J\f and J\f as 

/ \ 2 


»(B) = tmir = t 

1=1 1=1 7=1 


^ i¥i 


Let us now consider the function: 


/(v) = E 


tyivtyjv 


\h-hj 

i¥i 


(50) 


(51) 


Assuming D/v" diagonal and with constant diagonal elements k, we 
have 


k [ f(v)dv > d&(B ), 
Jb 


(52) 


in fact: 


a*(B) « kV V £ve3.M fry 4>/v 

* A; — A 7 


7=1 \ J=1 'V 

V 7V» 


* I E 


*. 7=1 hj 

J¥i 


= k £ /(v). (53) 

v£d¥A 
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